MEANWR1.doc.

' ANWR Results Program
'  Description: This is a Monte Carlo simulation program to
'    compute estimated in-place and recoverable oil and gas
'    resources by play.
'  Input: Input data is obtained from Excel 97 worksheets Oil,
'    and Play-Prospect.
'  Output: The principal output is specified in worksheet
'    Summary.  It consists of estimates of the mean and
'    standard deviation and selected quartiles of in-place
'    and recoverable oil, NGL, and gas.  Size-frequency
'    distributions are given in worksheet Distns.
'    Supplementary results used to check the program and data
'    are given in worksheet Supple
'    In addition, three ASCII files are created:
'     ProspData is detailed prospects data
'     PlayData is summary data by simulation run
'     HydroData is first 100 records-hydro output for testing
'  Operating instructions: In the Excel 97 Workbook go to
'    Tools, Macro, Macros.  Then run ANWR1_VBA.
'
' Written in Visual Basic for Applications (Excel 97)
'  John H. Schuenemeyer, USGS, 3/31/98-11:30 am

Option Explicit
'Declarations
Private user_sheet  As Object
Private data_sheet  As Worksheet
Private chart_range As Range
Private chart_sheet As Chart
Private outer_max   As Integer
Private outer_index As Integer
Private inner_max   As Integer

Sub ANWR1_VBA()
'Define variables
Dim i, j, m, ni, nj, fsc, os, minfs, njknt As Integer
Dim r12, orf, grf, fplay As Variant
Dim proboil, adg, adgr, fpros As Variant
Dim GOR, ag, gor1, gor2 As Variant
Dim gast, gastr As Variant
Dim so, sor, nglr1, nglr2 As Variant
Dim sngl, snglr, sg, sgr As Variant
Dim pd, fvf1, fvf2, fvf, depth, depthadd As Variant
Dim nglratio, ngl, nglr, nglt, ngltr As Variant
Dim nglrar As Variant, nglra As Variant
Dim nglrag As Variant, nglrag1 As Variant, nglrag2 As Variant
Dim temp, temp1, temp2 As Variant
Dim press, press1, press2, press3, press4 As Variant
Dim zc As Variant, zc1 As Variant, zc2 As Variant
Dim marray As Integer, seed As Integer
Dim OilR As Variant, OilA As Variant, GasR As Variant, GasA As Variant
Dim ProsDist As Variant, njr As Variant, TotClos As Variant
Dim Mclos As Variant, Ntotc As Integer, Pclos As Variant
Dim mdiv, playa, mplaya As Variant
Dim z(1 To 6), TotR(1 To 7), maxfsc(1 To 14)
Dim TotIP(1 To 7), NumDep(1 To 4), NDRep(1 To 2)
Dim SumIP(1 To 7) As Double, SumR(1 To 7) As Double
Dim Sum2IP(1 To 7) As Double, Sum2R(1 To 7) As Double
Dim tmpIP, tmpR As Variant, c(1 To 7)

'Obtain input values from worksheets

Worksheets("Play-Prospect").Activate
 Range("h18").Select
 fplay = Selection.Value  'prob of favorable play
 Range("h24").Select
 fpros = Selection.Value  'prob of favorable prospect
 Range("g30").Select
 proboil = Selection.Value 'prob prospect is oil
 Range("f3").Select
 minfs = Selection.Value   'cutoff value
 Range("j8:k14").Select
 ProsDist = Selection.Value 'num of prospects distn
 
Worksheets("Oil").Activate
  Range("c4").Select
   playa = Selection.Value  'Play area in thousand of acres
  Range("g19").Select
   depthadd = Selection.Value 'Surface to sea level
  Range("c22").Select         ' thousands of feet
   orf = Selection.Value    'oil rcovery factor, %
   orf = orf * 0.01
  Range("c26").Select
   fvf1 = Selection.Value   'fvf equation parameter
  Range("e26").Select
   fvf2 = Selection.Value   'fvf equation parameter
  Range("c29").Select
   gor1 = Selection.Value   'GOR equation parameter
  Range("e29").Select
   gor2 = Selection.Value   'GOR equation parameter
  Range("c31").Select
   nglr1 = Selection.Value  'ngl ratio eqn parameter
  Range("e31").Select
   nglr2 = Selection.Value  'ngl ratio eqn parameter
  Range("j10:u16").Select
   pd = Selection.Value     'hydro volume parameters
   
Worksheets("Gas").Activate
  Range("d19").Select
   grf = Selection.Value   'gas rcovery factor, %
   grf = grf * 0.01
  Range("c23").Select
   nglrag1 = Selection.Value   'ngl+cond to NA gas eqn parameter
  Range("e23").Select
   nglrag2 = Selection.Value   'ngl+cond to NA gas eqn parameter
  Range("d42").Select
   press1 = Selection.Value   'gas press eqn parameter
  Range("f42").Select
   press2 = Selection.Value   'gas press eqn parameter
  Range("d43").Select
   press3 = Selection.Value   'gas press eqn parameter
  Range("f43").Select
   press4 = Selection.Value   'gas press eqn parameter
  Range("d45").Select
   temp1 = Selection.Value   'gas temp eqn parameter
  Range("f45").Select
   temp2 = Selection.Value   'gas temp eqn parameter
  Range("d47").Select
   zc1 = Selection.Value   'gas compress eqn parameter
  Range("f47").Select
   zc2 = Selection.Value   'gas compress eqn parameter

Worksheets("Distns").Activate 'In place
  Range("b5:n6").Select
   so = Selection.Value     'store num of dep & total oil
  Range("b9:n12").Select    'gas
   sg = Selection.Value
  Range("b16:n18").Select   'ngl
   sngl = Selection.Value
  Range("b22:n23").Select   'Recoverable
   sor = Selection.Value    'store num of dep & in-place oil
  Range("b26:n29").Select
   sgr = Selection.Value    'gas-recoverable
  Range("b33:n35").Select   'ngl
   snglr = Selection.Value
Worksheets("Supple").Activate  'supplementary results
  Range("e5").Select
   Mclos = Selection.Value
   
'Clear arrays
marray = 13     'marray is the maximum number of log 2 classes
For j = 1 To marray
 For i = 1 To 2
  so(i, j) = 0
  sor(i, j) = 0
 Next i
 For i = 1 To 3
  sngl(i, j) = 0
  snglr(i, j) = 0
 Next i
 For i = 1 To 4
  sg(i, j) = 0
  sgr(i, j) = 0
 Next i
Next j
For j = 1 To 7
 SumIP(j) = 0
 SumR(j) = 0
 Sum2IP(j) = 0
 Sum2R(j) = 0
Next j
' The following array stores the maximum field size classes
' Order is oil,ngl,gas for in place and then recoverable
For i = 1 To 14
 maxfsc(i) = 0
Next i
For i = 1 To 4
 NumDep(i) = 0
Next i
Mclos = 0
Ntotc = 0

Open "ProspData" For Output As #1    ' Open file for output.
Write #1, "SimNum", "ProspNum", "O1G2", "O/G", "Closure", "Depth", "Other params"
Open "PlayData" For Output As #2
Write #2, "SimNum", "OilR", "AssDGas", "NA_Gas", "All_Gas", "NGL_ADG", "NGL_NAG", "All_NGL"
Open "HydroData" For Output As #3    ' Open file for output.
Write #3, "SimNum", "ProspNum", "O1G2", "O/G", "Thick", "Closure", "Porosity", "100-WS", "TrapFill"

m = 10000  'm is num of simulation runs=num successful plays
mdiv = Int(m / fplay + 0.5) 'total number of plays
seed = -374 'random number seed, negative integer
Rnd (seed)  'initializes uniform random number generator
' subtraction for computing variance if needed
c(1) = 0 '3500 'est for in place oil mean
c(2) = 0 '1900 'est for in place assoc diss gas mean
c(3) = 0 '10500 'est for in place NA gas mean
c(4) = c(2) + c(3) 'est for in place total gas mean
c(5) = 0 '1900 'est for in place ngl from assoc diss gas mean
c(6) = 0 '10500 'est for in place ngl from NA gas mean
c(7) = c(5) + c(6) 'est for in place total ngl mean

'main simulation loop
For i = 1 To m
 For j = 1 To 7
  TotR(j) = 0
  TotIP(j) = 0
 Next j
NDRep(1) = 0
NDRep(2) = 0
TotClos = 0
njknt = 0
Call SingleS(ProsDist, 1, njr)
 ni = CInt(njr)     'ni is number of prospects, CInt rounds
 
' Prospect loop
For nj = 1 To ni

' Check favorability of nj th prospect
If Rnd() < fpros Then
 njknt = njknt + 1   'count of successful prospects in play i
 ngl = 0
 nglra = 0
 nglr = 0
 nglrar = 0
 GasA = 0
 adg = 0
 GasR = 0
 adgr = 0
 
' Is reservoir oil or gas?
 If Rnd() < proboil Then
'OIL

 os = 1     'oil reservoir
  Call SingleS(pd, 11, depth)   'Sample depth
  depth = depth + depthadd
'In-place oil, mm bbl
'formation volume factor
  If depth < 2.17 Then     ' thousands of feet
   fvf = 1
  Else
   If depth > 12.15 Then   ' thousands of feet
    fvf = 1.5
     Else
    fvf = fvf1 + fvf2 * depth
   End If
  End If
  
  Call SampleOHVP(z, pd, fvf, minfs)  'get size >=cutoff
  OilA = z(6)
  Call SizeClass(OilA, maxfsc(1), marray, fsc)
  so(1, fsc) = so(1, fsc) + 1
  so(2, fsc) = so(2, fsc) + OilA
  NDRep(1) = NDRep(1) + 1
  TotIP(1) = TotIP(1) + OilA
  
'in place assoc-dissolved gas, billions of cu. ft.
  GOR = 10 ^ (gor1 + gor2 * depth)
  adg = OilA * GOR * 0.001  '.001 for billions of cu. ft.
  Call SizeClass(adg / 6, maxfsc(2), marray, fsc)
  sg(2, fsc) = sg(2, fsc) + adg
  TotIP(2) = TotIP(2) + adg
    
' in place ngl from assoc-diss gas, millions bbl
  If depth >= 15.674 Then
    nglratio = 100     'units=bbl/millions cu ft
    Else
    nglratio = 1000000 / (nglr1 * Exp(nglr2 * depth))
  End If
  ngl = adg * nglratio * 0.001 ' ngl, millions bbl
  Call SizeClass(ngl, maxfsc(5), marray, fsc)
  sngl(1, fsc) = sngl(1, fsc) + ngl
  TotIP(5) = TotIP(5) + ngl
  
 ' recoverable oil, mm bbl
 OilR = OilA * orf   ' orf recoverable oil
  Call SizeClass(OilR, maxfsc(8), marray, fsc)
  sor(1, fsc) = sor(1, fsc) + 1
  sor(2, fsc) = sor(2, fsc) + OilR
  TotR(1) = TotR(1) + OilR
  'wellpoil=.16*OilR*z(2) ???
    
' recoverable assoc-dissolved gas, billions of cu. ft.
  adgr = OilR * GOR * 0.001
  Call SizeClass(adgr / 6, maxfsc(9), marray, fsc)
  sgr(2, fsc) = sgr(2, fsc) + adgr
  TotR(2) = TotR(2) + adgr
  
' recoverable ngl from assoc-diss gas, millions bbl
  nglr = adgr * nglratio * 0.001
  Call SizeClass(nglr, maxfsc(12), marray, fsc)
  snglr(1, fsc) = snglr(1, fsc) + nglr
  TotR(5) = TotR(5) + nglr
 Else
 
 'GAS
 'Note that provision was made for sampling from
 ' separate distributions for non-associated gas
 ' however in this study all hydrocarbon distributions
 ' are specified on the Oil worksheet.
 os = 2     'code for gas
 Call SingleS(pd, 11, depth)    'Sample depth
 depth = depth + depthadd
 If press4 = 0 Or depth <= 10 Then
  press = press1 + press2 * depth 'Pressure, psi
  Else
  press = press3 + press4 * depth
 End If
 temp = temp1 + temp2 * depth
 Call gasc(press, temp, depth, zc)
 Call SampleGHVP(z, pd, zc, temp, press, minfs)
  GasA = z(6)
  
 ' in place NA gas, billions of cu. ft.
  Call SizeClass(GasA / 6, maxfsc(3), marray, fsc) 'size class
  sg(1, fsc) = sg(1, fsc) + 1
  sg(3, fsc) = sg(3, fsc) + GasA
  NDRep(2) = NDRep(2) + 1
  TotIP(3) = TotIP(3) + GasA
  
'  in place natural gas liquids plus condensate to
'    non-associated gas (bbls/million cf)
  If depth >= 15.674 Then
    nglratio = 100       'units=bbl/million cf
    Else
    nglratio = 1000000 / (nglr1 * Exp(nglr2 * depth))
  End If
  nglra = GasA * nglratio * 0.001  'ngl in mm bbl
  Call SizeClass(nglra, maxfsc(6), marray, fsc) 'size class
  sngl(2, fsc) = sngl(2, fsc) + nglra  'ngl from NA gas
  TotIP(6) = TotIP(6) + nglra 'contribution to ngl from NA gas
  
 ' recoverable NA gas
  GasR = GasA * grf   'recoverable gas, billions cu. ft.
  Call SizeClass(GasR / 6, maxfsc(10), marray, fsc)
  sgr(1, fsc) = sgr(1, fsc) + 1
  sgr(3, fsc) = sgr(3, fsc) + GasR
  TotR(3) = TotR(3) + GasR
  'wellprodg=.62*GasR*z(2)
  
  'recov natural gas liquids plus condensate to
  '  non-associated gas(bbls/million cf)
  nglrar = GasR * nglratio * 0.001  'ngl in mm bbl
  Call SizeClass(nglrar, maxfsc(13), marray, fsc) 'size class
  snglr(2, fsc) = snglr(2, fsc) + nglrar  'ngl from NA gas
  TotR(6) = TotR(6) + nglrar 'contribution to ngl from NA gas
  
 End If     'end of oil/gas computations
 
 'get total gas and total ngl
 gast = adg + GasA      'total in-place gas
 If gast > 0 Then
  Call SizeClass(gast / 6, maxfsc(4), marray, fsc) 'size class
  sg(4, fsc) = sg(4, fsc) + gast
  TotIP(4) = TotIP(4) + gast
  gastr = adgr + GasR    'total recoverable gas
  Call SizeClass(gastr / 6, maxfsc(11), marray, fsc) 'size class
  sgr(4, fsc) = sgr(4, fsc) + gastr
  TotR(4) = TotR(4) + gastr
 End If
 nglt = ngl + nglra     'total in-place ngl
 If nglt > 0 Then
  Call SizeClass(nglt, maxfsc(7), marray, fsc) 'size class
  sngl(3, fsc) = sngl(3, fsc) + nglt
  TotIP(7) = TotIP(7) + nglt
  ngltr = nglr + nglrar  'total recoverable ngl
  Call SizeClass(ngltr, maxfsc(14), marray, fsc) 'size class
  snglr(3, fsc) = snglr(3, fsc) + ngltr
  TotR(7) = TotR(7) + ngltr
 End If
 TotClos = TotClos + z(2)
 If i <= 100 Then
  If os = 1 Then
   Write #3, i, njknt, os, z(6), z(1), z(2), z(3), z(4), z(5)
  Else
   Write #3, i, njknt, os, z(6), z(1), z(2), z(3), z(4), z(5)
  End If
 End If
 If os = 1 Then
  Write #1, i, njknt, os, OilR, z(2), depth, adgr, nglr, fvf, GOR, nglratio
 Else
  Write #1, i, njknt, os, GasR, z(2), depth, ngltr, press, temp, zc, nglratio
 End If
End If
' no reservoir
Next nj   'end of prospect loop

If TotClos > Mclos Then Mclos = TotClos
If TotClos > playa Then Ntotc = Ntotc + 1

Write #2, i, TotR(1), TotR(2), TotR(3), TotR(4), TotR(5), TotR(6), TotR(7), TotIP(1), TotIP(2), TotIP(3), TotIP(4), TotIP(5), TotIP(6), TotIP(7)
For j = 1 To 7
 tmpIP = TotIP(j) - c(j)
 tmpR = TotR(j) - c(j) * 0.5
 SumIP(j) = SumIP(j) + tmpIP
 SumR(j) = SumR(j) + tmpR
 Sum2IP(j) = Sum2IP(j) + tmpIP * tmpIP
 Sum2R(j) = Sum2R(j) + tmpR * tmpR
Next j
For j = 1 To 2
 NumDep(j) = NumDep(j) + NDRep(j)
 NumDep(j + 2) = NumDep(j + 2) + NDRep(j) * NDRep(j)
Next j

Next i    'end of simulation loop

'Summary statistics
' get averages by size class

For i = 1 To 7
 Sum2IP(i) = Sqr((Sum2IP(i) - (SumIP(i) / mdiv) * SumIP(i)) / (mdiv - 1))
 Sum2R(i) = Sqr((Sum2R(i) - (SumR(i) / mdiv) * SumR(i)) / (mdiv - 1))
 SumIP(i) = SumIP(i) / mdiv + c(i)
 SumR(i) = SumR(i) / mdiv + c(i) * 0.5
Next i
For j = 1 To 2
 NumDep(j + 2) = Sqr((mdiv * NumDep(j + 2) - NumDep(j) * NumDep(j)) / ((mdiv - 1) * mdiv))
 NumDep(j) = NumDep(j) / mdiv
Next j
For j = 1 To marray
 For i = 1 To 2
  so(i, j) = so(i, j) / mdiv
  sor(i, j) = sor(i, j) / mdiv
 Next i
 For i = 1 To 3
  sngl(i, j) = sngl(i, j) / mdiv
  snglr(i, j) = snglr(i, j) / mdiv
 Next i
 For i = 1 To 4
  sg(i, j) = sg(i, j) / mdiv
  sgr(i, j) = sgr(i, j) / mdiv
 Next i
Next j

Pclos = 100 * Ntotc / m

'Output statistics
Worksheets("Summary").Activate 'in place mean and std dev
  Range("b5").Select
   Selection.Value = SumIP(1)
  Range("c5").Select
   Selection.Value = Sum2IP(1)
  Range("b7").Select
   Selection.Value = SumIP(2)
  Range("c7").Select
   Selection.Value = Sum2IP(2)
  Range("b8").Select
   Selection.Value = SumIP(3)
  Range("c8").Select
   Selection.Value = Sum2IP(3)
  Range("b9").Select
   Selection.Value = SumIP(4)
  Range("c9").Select
   Selection.Value = Sum2IP(4)
  Range("b11").Select
   Selection.Value = SumIP(5)
  Range("c11").Select
   Selection.Value = Sum2IP(5)
  Range("b12").Select
   Selection.Value = SumIP(6)
  Range("c12").Select
   Selection.Value = Sum2IP(6)
  Range("b13").Select
   Selection.Value = SumIP(7)
  Range("c13").Select
   Selection.Value = Sum2IP(7)
  Range("b17").Select            'recoverable mean and std dev
   Selection.Value = SumR(1)
  Range("c17").Select
   Selection.Value = Sum2R(1)
  Range("b19").Select
   Selection.Value = SumR(2)
  Range("c19").Select
   Selection.Value = Sum2R(2)
  Range("b20").Select
   Selection.Value = SumR(3)
  Range("c20").Select
   Selection.Value = Sum2R(3)
  Range("b21").Select
   Selection.Value = SumR(4)    ' Mean NA Gas
  Range("c21").Select
   Selection.Value = Sum2R(4)
  Range("b23").Select
   Selection.Value = SumR(5)
  Range("c23").Select
   Selection.Value = Sum2R(5)
  Range("b24").Select
   Selection.Value = SumR(6)    ' Mean NA Gas
  Range("c24").Select
   Selection.Value = Sum2R(6)
  Range("b25").Select
   Selection.Value = SumR(7)
  Range("c25").Select
   Selection.Value = Sum2R(7)
  Range("b28").Select
   Selection.Value = NumDep(1)  ' Number of pools, mean and std.dev
  Range("c28").Select
   Selection.Value = NumDep(3)
  Range("b29").Select
   Selection.Value = NumDep(2)
  Range("c29").Select
   Selection.Value = NumDep(4)

  Range("b32").Select
   Selection.Value = m      'number of simulations
  Range("b33").Select
   Selection.Value = mdiv   'number of play runs
   
Worksheets("Distns").Activate
  Range("b5:n6").Select     'in place
   Selection.Value = so     'num and oil
  Range("b9:n12").Select
   Selection.Value = sg     'num NA gas, ass dis gas & NA gas
  Range("b16:n18").Select
   Selection.Value = sngl   'ngl from ass dis gas & NA gas
  Range("b22:n23").Select   'recoverable
   Selection.Value = sor    'num and oil
  Range("b26:n29").Select
   Selection.Value = sgr    'num NA gas, ass dis gas & NA gas
  Range("b33:n35").Select
   Selection.Value = snglr  'ngl from ass dis gas & NA gas
    
Worksheets("Supple").Activate
  Range("e5").Select
   Selection.Value = Mclos   'Max area of closures
  Range("e6").Select
   Selection.Value = Pclos   '%total closure>play area|success
  Range("e13:r13").Select
   Selection.Value = maxfsc  'Max field size classes
  Range("e8").Select
   Selection.Value = marray  'Max num log 2 size classes
  Range("e3").Select
   Selection.Value = seed    'Random number seed
Close #1
Close #2
Close #3

End Sub
Sub SampleOHVP(z, pd, fvf, minfs)
' Computes Oil in place >= cutoff
'  Note correlation between z(3) [porosity] and
'   z(4)[1-ws] is one.
Dim prod As Variant, j As Integer, k As Integer, u As Variant
Dim j2 As Integer
Line1:
 prod = 1
  For j = 1 To 5
    If j <> 4 Then u = Rnd
    j2 = j * 2 - 1
     For k = 2 To 7
       If u > pd(k - 1, j2 + 1) And u <= pd(k, j2 + 1) Then
         z(j) = pd(k - 1, j2) + ((u - pd(k - 1, j2 + 1)) / (pd(k, j2 + 1) - pd(k - 1, j2 + 1))) * (pd(k, j2) - pd(k - 1, j2))
         prod = prod * z(j)
       End If
     Next k
  Next j
' in place oil in millions of bbl
z(6) = prod * 7.758 * 10 ^ -6 / fvf
If z(6) >= minfs Then
 Exit Sub
Else
 GoTo Line1
End If
End Sub
Sub SampleGHVP(z, pd, zc, temp, press, minfs)
' Computes Gas in place >= cutoff
'  Note correlation between z(3) and z(4) is one.
Dim prod As Variant, j As Integer, k As Integer, u As Variant
Dim j2 As Integer
Line1:
 prod = 1
  For j = 1 To 5
    If j <> 4 Then u = Rnd
    j2 = j * 2 - 1
     For k = 2 To 7
       If u > pd(k - 1, j2 + 1) And u <= pd(k, j2 + 1) Then
         z(j) = pd(k - 1, j2) + ((u - pd(k - 1, j2 + 1)) / (pd(k, j2 + 1) - pd(k - 1, j2 + 1))) * (pd(k, j2) - pd(k - 1, j2))
         prod = prod * z(j)
       End If
     Next k
  Next j
'in place gas in billions of cu ft
z(6) = prod * 0.00000154 * (press / temp) / zc
' 1540e-09=43.56 * 10 ^ -9 * (519.7 / 14.7)
If z(6) / 6 >= minfs Then
 Exit Sub
Else
 GoTo Line1
End If
End Sub
Sub SingleS(pd, j, z)
' Random Deviate from Fractile Distribution
Dim k As Integer, u As Variant

   u = Rnd
    For k = 2 To 7
       If u > pd(k - 1, j + 1) And u <= pd(k, j + 1) Then
         z = pd(k - 1, j) + ((u - pd(k - 1, j + 1)) / (pd(k, j + 1) - pd(k - 1, j + 1))) * (pd(k, j) - pd(k - 1, j))
       End If
    Next k
End Sub
Sub SizeClass(size, maxf, marray, fsc)
' Computes size class in log based 2
 fsc = Int(Log(size) / Log(2)) - 1 ' size class
  If fsc < 1 Then fsc = 1
  If fsc > maxf Then maxf = fsc
  If fsc > marray Then fsc = marray
 End Sub
Sub gasc(p0, t0, depth, z)
'Estimates gas compressibility
'  Ref: Microcomputer Programs for Petroleum Engineers
'  p- pressure
'  t- degrees Rankine
'  depth- 1000 ft
'  z- gas compressibility

Dim dr(1 To 6), drp(1 To 6)
Dim i, k As Integer
Dim a, b, c, e, f, g, t2, fdr, gdr, t, p As Variant
If depth <= 3 Then      'depth <=3,000 ft
  z = 1 - 0.11 * depth + 0.0125 * depth * depth
  Exit Sub
End If
' Use Mricocompute algorithm for depths > 3000 ft
' t and p are reducted temperature and pressure
'  using gas gravity of 0.6 and Fig 17 furnished by
'  John Quinn
t = t0 / 358
p = p0 / 672
a = 0.06423
f = 0.6845
dr(1) = 0.27 * p / t

b = 0.5353 * t - 0.6123
t2 = t * t
c = 0.3151 * t - 1.0467 - (0.5783 / t2)
e = 0.6816 / t2
g = 0.27 * p
For i = 1 To 5
drp(1) = dr(i)
For k = 2 To 6
 drp(k) = drp(k - 1) * drp(1)
Next k
fdr = a * drp(6) + b * drp(3) + c * drp(2) + t * drp(1) + e * drp(3) * (1 + f * drp(2)) * Exp(-f * drp(2)) - g
gdr = 6 * a * drp(5) + 3 * drp(2) + 2 * c * drp(1) + t + e * drp(2) * (3 + f * drp(2) * (3 - 2 * f * drp(2))) * Exp(-f * drp(2))
dr(i + 1) = dr(i) - fdr / gdr
Next i
z = dr(6)
End Sub
1


1